/**
  ETFDAQ Project -- online from ETFAna Online by SUN Yazhou
  \class TAMath
  \brief Math class, to provide some general math methods.
  \author SUN Yazhou, aisa.rabbit@163.com
  \date Created: 2020/07/09
  \date Last modified: 2023/10/18 by SUN Yazhou
  \copyright 2020-2023 SUN Yazhou
  \copyright Copyright (c) 2021-2024 Anyang Normal U. && IMP-CAS with LGPLv3 LICENSE
  \copyright MOMD project, Anyang Normal University, IMP-CAS
*/

#ifndef TAMath_h
#define TAMath_h

#include <cmath>
#include <string>
#include <iostream>

using std::string;
using std::cout;
using std::endl;

class TAMath{
public:
  TAMath(){}
  virtual ~TAMath(){}

  // physical functions
  static constexpr double Pi(){ return 3.14159265358979323846; }
  /// golden cut ratio
	static constexpr double Alpha(){ return 0.61803398874989484820458683436565; }
  /// rad per degree
	static constexpr double DEGREE(){ return 0.01745329251994329547; }
	static constexpr double Sqrt2(){ return 1.414213562373095048802; }
	static constexpr double Sqrt3(){ return 1.73205080756887719318; }
  static constexpr double uMeV(){ return 931.494061; } ///< u (C12=12.u) in MeV
  static constexpr double hbarc(){ return 197.32705; } ///< \retval hbar*c in MeV*fm
  static constexpr double FineStructureConstant(){ return 7.29735307964482e-3; } ///< 1./137.0359895
  static constexpr double Pion0Mass(){ return 134.9770; } ///< \retval pion0 mass in MeV
  static constexpr double PionPMMass(){ return 139.5706; } ///< \retval pion+- mass in MeV
  static constexpr double KaonPMMass(){ return 493.677; } ///< \retval pion+- mass in MeV
  static constexpr double c0(){ return 299.792458; } ///< \reteval light speed in mm/ns
  // CC = 0.321840605 = e0/(u0*c0*1E6) SI unit
  static constexpr double CC(){ return 0.32184060506023993; }
  static constexpr double MB(){ return 1048576.; } // = 1024*1024

  /// \note r_global = (R.r_local)+r0
	static void rotate(double *pOut, const double *pIn, const double *angIn);
  /// distance of a point to a 2-D straight line
  /// l: x = k*z+b or k*z-x+b=0
  static double DistanceToLine(double z, double x, double k, double b){
    return (k*z-x+b)/sqrt(1.+k*k);
  } // end member function DistanceToLine
  /// distance between 2 points in zOx plane
  static double Distance(double z0, double x0, double z1, double x1){
    return sqrt(sqr(z0-z1) + sqr(x0-x1));
  } // end member function DistanceToLine

  /// from Ek per u to velocity for a generic nucleon of mass uMeV
  static double EkPerUToBeta(double pEk){
    const double r = 1. + uMeV()/pEk;
    return sqrt(2.*r-1.)/r;
  } // end member function EkPerUToBeta
  /// from Ek per u to momentum for a generic nucleon of mass uMeV
  static double EkPerUToP(double pEk){
    const double b = EkPerUToBeta(pEk);
    return b/sqrt(1.-b*b)*uMeV();
  }
  static double WoodsSaxon(double r, double a0, double r0){ return 1./(1.+exp((r-r0)/a0)); }

  /// sign of a number
  static double sign(double c){
    if(c >= 0.) return 1.;
    return -1.;
  }
  /// sign of a number
  static double sign(int c){ return sign(double(c)); }
  /// \retval (-)^m
  static double minus(int n){ return n & 1 ? -1 : 1; }
  static double sqr(double a){ return 0. == a ? a : a*a; }
  static double innerProd(const double *p0, const double *p1, int n = 3){
    double prod = 0.;
    for(int i = n; i--;) prod += p0[i]*p1[i];
    return prod;
  } // end member function innerProd
  static double norm(const double *p, int n = 3){ return sqrt(innerProd(p, p, n)); }
  /// \retval computes (a^2+b^2)^0.5 without desctructive underflow or overflow
  static double norm(double a, double b){
    a = fabs(a); b = fabs(b);
    if(a > b) return a*sqrt(1.+sqr(b/a));
    else return 0. == b ? 0. : b*sqrt(1.+sqr(a/b));
  } // end member function norm
  static double infNorm(int n, const double *f){
    double m = 0.;
    for(int i = n; i--;) if(fabs(f[i]) > m) m = fabs(f[i]);
    return m;
  } // end member function InfNorm
  // generic programming :)
  /// \retval: sum of square of each one in the parameter list
  static double Sum2(){ return 0.; }
  template<class T, class ...Args> // typename preferred, class deprecated
  static T Sum2(const T &v, const Args &...args){
    return v*v + Sum2(args...);
  }
  template<class T>
  static bool Within(const T &v, const T &A, const T &B){
    if(v >= A && v <= B) return true;
    if(v >= B && v <= A) return true;
    return false;
  }
  static double ang(double k){
    return atan(k) / DEGREE();
  } // end member function

  ///< print array in an orderly fashion
  template<class T>
  static void peek(const T *v, int len){
    for(int i = 0; i < len; i++){
      // cout << v[i] << '\t';
      cout << int(v[i]) << '\t'; // DEBUG
      if(i != 0 && i % 10 == 0) cout << endl;
    } // end for
    cout << "\033[32;1m---\033[0m" << endl;
  } // end member function peek

  static bool endwith(const string &s, const string &ss){
    if(s.size() < ss.size()) return false;
    return s.substr(s.size()-ss.size(), ss.size()) == ss;
  } // end member function endwith
  static bool beginwith(const string &s, const string &ss){
    if(s.size() < ss.size()) return false;
    return s.substr(0, ss.size()) == ss;
  } // end member function beginwith
  static bool has(const string &s, const char *p){
    return s.find(p) != string::npos;
  } // end member function has
  // strcmp, but ignore case
  static bool istrcmp(const string &s0, const string &s1);
  static void uppercase(string &s);

  static long filesize(const string &f);
  static double filesizeInMB(const string &f);

  /// linear interpolation
  /// \param kk stores the slope upon request
  static double lii(double x0, double y0, double x1, double y1, double x,
    double *kk = nullptr);
  static double lii(const double *p0, const double *p1, double x, double *k = nullptr){
    return lii(p0[0], p0[1], p1[0], p1[1], x, k);
  } // end member function lii

  /// \retval calculate combination number
  /// returns the value ln[Gamma(xx)] for xx > 0.
  static double gammaln(double xx);
  // returns ln(n!)
  static double factln(int n);
  /// \retval n!
  static double Factorial(int n);
  /// \param n, m: n should not be smaller than m
  static int Binomial(int n, int m);
  /// beta function, Binomial coefficients with float arguments
  static double Beta(double z, double w);
  /// \retval n!!
  static int BiFactorial(int n);
  /// the incomplete gamma function gammap, gammpa(a,x)=P(a,x)=gamma(a,x)/Gamma(a)
  static double gammap(double a, double x);
  /// the incomplete gamma function gammap, gammpa(a,x)=P(a,x)=gamma(a,x)/Gamma(a)
  static double gammaq(double a, double x);
  /// returns the incomplete gamma function P(a,x) evaluated by its series
  /// representation. Optionally returns ln[Gamma(a)] as gln.
  static double gser(double a, double x, double *gln = nullptr);
  /// returns the incomplete gamma function P(a,x) evaluated by its continued fraction
  /// representation. Optionally returns ln[Gamma(a)] as gln.
  static double gcf(double a, double x, double *gln = nullptr);

  static double BetaGamma(double beta);

  static void normalize(double *p, int n = 3);
};

#endif
